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ABSTRACT 

Aims. We explore the effects of photodissociation at the stages of post-asymptotic giant branch stars to find a mechanism able to 
produce multi-polar shapes. 

Methods. We perform two-dimensional gasdynamical simulations to model the effects of photodissociation in proto-planetary nebu- 
lae. 

Results. We find that post-asymptotic giant branch stars with ~ 7, 000 K or hotter are able to photodissociate a large amount of 
the circumstellar gas. We compute several solutions for nebulae with low-velocity multi-lobes. We find that the early expansion of a 
dissociation front is crucial to understand the number of lobes in proto-planetary nebulae. 

Conclusions. A dynamical instability appears when cooling is included in the swept-up molecular shell. This instability is similar 
to the one found in photoionization fronts, and it is associated with the thin-shell Vishniac instability. The dissociation front exac- 
erbates the growth of the thin-shell instability, creating a fast fragmentation in shells expanding into media with power-law density 
distributions such as r~ 2 . 
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1. Introduction 

Several surveys of proto-planetary nebulae (PPNs) and plan- 
etary nebulae (PNs) (|Sahai et al. (2007)1 |Chu et al. (1987) 
Balick (1987)1 |Schwarz et al. (1992)1 |Stanghellini et al. (1993) 



Manchado et al. (1996a)[ |Manchado et al. (1996b)[ ) show a 
rich variety of shapes, and they have been cataloged in a 
series of morphological classes: bipolar, elliptical, point- 
symmetric, irregular, spherical, quadrupolar, and multi-polar. 
While the first five former classes have been reproduced 
with different scenarios by hydrodynamical and magnetohy 
drodynamical simul ations \ Icke (1988)| 



Soker & Livio (1989), Melle ma et al. (1991)1 |Icke et al. (1992) 
Frank & Mellema (1994)1 |Dwarkadas et al. (19967 
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Ickeet al. (1989) 



Garcfa-Segura (1997) 



Rijkhorst et al. (2005) 



Rozyczka & Franco (1996) 
Garcfa-Segura et al. (1999l{ 
Garcfa-Segura et al. (2005)) , the last two former classes 
have eluded any attempt of modeling, except the one by 
|Rijkhorst et al. (2005)| and |Matt et al. (2006)| In general, all 
previous scenarios are based on the collimation of winds. 
Two examples of low-expanding, multi-polar PPNs are IRAS 
19024+0044 (ISahai et al. (2007b)) and IRAS 19475+3119 
(Sahaietal. (2007c)J). It is interesting that these two nebulae 
neither show any kinematic fingerprint in their optical spectra 
nor a perfect symmetry. 

Only recently, the hypothesis that bipolar and multipolar 
lobes of PNs are the result of ionization and illumination was 
proposed by |Kwok (2010)[ He suggested that multipolar nebu- 
lae with similar lobe sizes are not caused by simultaneous ejec- 
tion of matter in several directions, but by the leakage of UV 
photons into those directions. Therefore, the explanation to the 
multipolar phenomenon does not lie in the dynamical ejection, 
but in how such holes were formed. The similarity in sizes of the 



multiple lobes was therefore naturally explained, and no simul- 
taneous ejection was required. 

According to Kwok's idea, it is necessary to find a 
mechanism that is able to create holes in the surrounding 
medium of PPNs, different to stellar winds or collimated out- 
flows of any kind, such as jets originated by accretion disks 
(Livio & Pringle (1997)), or the effects from photoionization 
(Garcfa-Segura & Franco (1996)), because post-asymptotic gi- 
ant branch stars (post- AGB s) are still not hot enough to produce 
a considerable ionization of the medium. 

| Diaz-Miller, Franco & Shore (1998)] made an important 
contribution to the study of photodissociation regions (PDRs) 
around main-sequence stars. They outlined the importance of 
PDRs for low-temperature stars (down to 7,500 K), because the 
number of non-ionizing photons in the range of 912 A< A < 
1100 A is several orders of magnitude larger than the ionizing 
photons ( A < 912 A), which are able to dissociate the outlying 
molecular gas, forming PDRs. 

A number of previous studies ( |Latter et al. (2000) 
and references therein) have outlined the importance of 
PDRs in the context of PNs. However, previous models 
(Natta & Hollenbach (1998)) have focused on stellar tempera- 
tures above 30,000 K , or even higher, such as the case of the 
Helix Nebula ( |Henney et al. (2007)) . 

Inspired by Kwok's idea, we explore the effects of pho- 
todissociation at the stages of post-AGBs, when stars still 
have low temperatures (~ 7,000 K), i.e., they are still in 
the so-called transition time (5,000 < 10,000 K), and it is 
found that PDRs are able to explain several observational fea- 
tures found in PPNs. The transition time of post-AGBs is not 
yet fully understood ( Stanghellini & Renzini (2000)), since the 
mass-lost rate at these stages is basically unknown . Stellar 
evolution calculations of PN central stars usually avoid this 
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Fig. 1. Snapshots of gas density at 1, 2, 3, 4 (top), and 5, 6, 7, 8 xlO 3 yr (bottom) (log(min)=-22, log(max)=-18). 



stage and start at 10,000 K (Vassiliadis & Wood (1994)) , al- 
though there is a model for 3za msM q that starts at ~ 6,000 K 
(Bloc ker & Schonberner (1990)) . In the context of the present 
paper, long transition times > 10,000 yr (low mass stars) will 
produce dynamically, developed PDRs, while transition times 
shorter than 1,000 yr (massive stars) do not give enough time 
to develop a structured PDR and will form rapidly photoionized 
nebulae. 

In this paper, we use the study by 
| Diaz-Miller, Franco & Shore (1998)] as a starting point, 
and merely replace the constant density condition by a nu- 
merical integration for any density distribution. The paper is 
structured as follows: Sect. 2 describes the numerical models 
and physical approaches used in this study. Sect. 3 presents 
two-dimensional results of PPNs. Sect. 4 discusses the results 
and gives a summary. 



2. Numerical models and physical approaches 

The simulations have been performed with the magneto- 
hydrodynamical code ZEUS-3D version 3.4 (Clarke (1996)), 
developed by M. L. Norman and the Laboratory for 
Computational Astrophysics. This is a finite-difference, fully 
explicit Eulerian code descended from the code described 
in |Stone & Norman (1992)[ ZEUS-3D does not include radi- 
ation transfer, but we have implemented a simple approxi- 
mation ( | Diaz-Miller, Franco & Shore (1998)) to derive the lo- 
cation of the dissociation front for arbitrary density distribu- 
tions, similar to the approximation used to locate the ion- 
ization front (Tenorio-Tagle (1979) Bodenhei mer et al. (1979)| 
|Franco etal.(1989)| 



Garcfa-Segura & Franco (1996) i. This is 



done by assuming that dissociation equilibrium holds at all 
times, and that the gas is fully dissociated inside the PDR. 
One can define the PDR equilibrium radius similarly to the 
Stromgren radius ( j Diaz-Miller, Franco & Shore (19 98)), and 
the position of the dissociation front in any given direction {8, <f>) 
from the photodissociating source is given by 



An 



n lot — 

1.37 x 
0.15 



where n H o is the density of neutral hydrogen, 
n H ° + 2«h, is the total proton density, aj = 
1Q-' 7 (|Hollenbach&Salpeter(1971)|l, < p > 
( [Abgrali et al. (1992)| |Draine & Bertoldi (1996)) , and S D the 
rate of dissociating photons. In the context of this paper, we as- 
sume that Run = 0. The left side of Eq. (1) is the number of H2 
molecules formed in a given volume, while the right side is the 
total number of photodissociations per unit time, where < p > 
is the average dissociation probability for the wavelength range 
912-1100 A. 

The stellar parameters are taken from a model for Mzams = 
2.5M (Vassiliadis & Wood (1994)), and using a black-body ap- 
proximation similar to Villaver et al. (2002) for T e g = 6 - 7 — 8 - 
10 x 10 3 K, we obtain log(S D ) ~ 42,43,44,45s-'. 

The models include the cooling curves by 
Dalgarno & McCray (1972)] and |MacDonald & Bailey (198T) 



The cooling cutoff temperature of the molecular gas is set to 10 
K. Finally, the photodissociated gas is always kept at 10 3 K , so 
no cooling curve is applied to the PDRs (unless there is a shock 
inside the PDR). 

Equatorial density enhancements are introduced 
in the two-dimensional models using Eq. 2 to 5 in 

based on the wind-compression 
equations by Bjorkman & Cassinelli (1993) Here , Q. = v rot /v cr i t 
gives the proximity of the star to critical rotation (see Eq. 6 



«h° n wt a/r 2 dr — < p > So, 



(1) 



in |Garcfa-Segura et al. (1999) 1. We assume that critical ro- 
tation (or close) is achieved at some point in the blueward 
excursion (Lange r~(1998)| |Heger & Langer (1998)) during 
the transition time, similar to the model computed for 12 
M Q (Chita et al. (2008)). Common envelope scenarios probably 
produce similar effects, although there are not yet any available 
equations in the literature. 



3. Two-dimensional results 

The first two-dimensional model consists of three stages. In 
the initial stage, we initialized the computational domain with 
a spherical AGB wind with M = 1 x 10~ 5 M o yr 1 and 
Vco - lOkms -1 . Our grid consists of 200 x 180 equidistant 
zones in r and 0, respectively. The innermost radial zone lies at 
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Table 1. Sample of models 



0.10 



model 


M agb (Me yr 1 ) 


n 


80-a 


0.60 x 10- 5 


0.80 


80-b 


0.50 x 10~ 5 


0.80 


50-a 


0.40 x 10" 5 


0.50 


50-b 


0.38 x 10~ 5 


0.50 



0.05 



0.00 



r; = 2.725 x 10~ 3 pc , and the outermost zone at r Q = 0.109 pc . 
The angular extent is 90°. 

In the second stage, we introduce a second wind at the inner 
boundary with the same parameters, but it is assumed that the 
star is near critical rotation with O = 0.98. This is done for a 
period of 10 4 yr. We do not include here a smooth transition 
(from Q = — > 0.98) for simplicity. These two stages are then 
used as initial conditions for a third stage. We set all velocities in 
the grid to zero. This is done to avoid the expansion of the AGB 
wind out of the grid and to compute only the expansion of the 
PDR. 

The first snapshot in Fig. 1 (at 1,000 yr) shows the structure 
generated by both AGB winds, an outer spherical wind almost 
at the edge of the grid, followed by a bipolar wind. This bipo- 
larity is owing to the smaller values of the radial velocity toward 
the equator in the AGB wind with critical rotation (see Eq. 2 in 
Garcf a-Segura et al. (1999)} , which produces an apparent empty 



space because of a rarefaction wave in the wind (dark areas in the 
plot). The dark areas would be filled with gas using a transition, 
evolving wind with Q. — — > 0.98. 

In the third stage, we simulate the conditions of a post- AGB 
wind near critical rotation (Q = 0.98), decreasing the mass-loss 
rate down to M = 1.73 x 1O~ 7 M yr~' . Here, we assume that 
the star has reached 7,000 K , and the corresponding UV pho- 
tons reach So — 10 43 s -1 . As before, we artificially set the wind 
velocity to zero, just to compute only the contribution from the 
thermal expansion of the photodissociated gas. By doing this, we 
are also assuming at the same time that the post- AGB star is not 
able to drive an efficient wind, because the stellar temperature is 
not high enough to produce a line-driven wind. 

Figure 1 shows the evolution of the gas density during 8,000 
yr. During the first thousand years, the lower density at the polar 
region allows the PDR to expand in the form of beams, digging 
low-density tunnels, preferentially at the polar axis (first snap- 
shot in Fig. 1). The gas is evacuated and pushed sideways in 
those beams because of the sudden increase of the temperature 
(from 10 to 1,000 K) in the PDR. During the second thousand 
years (second snapshot in Fig. 1), the growth of the beam at the 
polar axis ends because of the lateral expansion of the two ad- 
jacent beams, blocking the UV field owing to the increase of 
density at the axis. At the end of this phase at 2,000 yr, two 
thin lobes (0.08 pc long) are formed at each hemisphere (second 
snapshot in Fig. 1), producing a quadrupolar shape. 

Meanwhile, the density gradient between pole and equator in 
the PDR close to the central star produces a flow directed toward 
the polar axis. This bipolar flow converges at high latitudes, and 
is able to increase the density up to the point of producing a self- 
blocking of the UV field between 2,000 and 3,000 yr. At this 
moment, the PDR contracts and becomes confined by this flow 
at ~ 0.01 pc, and the external lobes cool down to form hydrogen 
molecules. After 3,000 yr, the PDR again expands in the polar 
direction in a bipolar or quadrupolar form (the dissociation front 
has funnel shape). 

The last snapshot in Fig. 1 (at 8,000 yr) shows the final devel- 
opment of our model. At 0.02 pc, a high-density blob produced 
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Fig. 2. From left to right, top to bottom: gas density of models 
80-a, 80-b, 50-a, 50-b at 9,000 yr of the evolution (log(min)=-22, 
log(max)=-18). Model numbers relate the approach to critical 
rotation (in percent). 



by the convergence of two oblique shocks defines the position 
of the dissociation front at the polar axis. Above this blob, the 
remnant lobes collimate the expansion of the molecular flow be- 
yond the dissociation front. At 0.045 pc, the molecular outflow 
reaches values of 8.7 km s _1 . 

The initial expansion of a dissociation front could be very 
sensitive to the conditions found at the density ramp. To illus- 
trate this behavior, we have computed a sample of four models 
in which we have changed the equator to pole density ratio by 
changing the stellar rotation (Q = 0.80, 0.50 ) , with a slightly 
different mass-loss rate for those with the same rotation (Table 
1). The other computational inputs have the same properties as 
the former model, except that we skip the first stage (the initial 
spherical wind), and directly initialize the computational mesh 
with the wind near critical rotation. The post-AGB wind is not 
included either, only the UV field. The results are shown in Fig. 
2, with snapshots at 9,000 yr. The angular extent is now 180° 
in 6, with 360 zones. We have introduced a 1% random noise 
in the initial density to allow for differences between the north- 
ern and the southern hemispheres, to emphasize that the lobes 
could reach different sizes and be formed at different angles, i.e., 
they are not necessarily symmetric. The expansion velocity of 
the lobes in Fig. 2 is lower than 1 km s _1 with respect to the 
AGB wind (= in the simulations). 
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We cannot discuss each model in detail here, but give a short 
summary. We find that the early expansion of a dissociation front 
is crucial to understand the number of lobes in the models, and 
we also find that the number of lobes is inversely proportional to 
the equator/pole density ratio. A PDR is able to expand in more 
directions when the ratio is low. For example, in the extreme 
case of Fig. 1, the PDR only occupies a few cells in the polar 
region, and the dissociation front digs a tunnel only through the 
polar axis (small solid angle), while in the case of model 50- 
b (Fig. 2), the initial PDR occupies many cells down to much 
lower latitudes, and the PDR expands into a larger solid angle. 
In the ideal case of a spherical wind (not shown in the figure), the 
two-dimensional result has the shape of a daisy flower, because 
the PDR expands in all directions. Note that the development of 
instabilities in a dissociation front is crucial to understand the 
expansion of PDRs (see Sect. 4). 

4. Discussion and summary 



Garcfa-Segura & Franco (1996) noted that a D-type dissoci 



ation front preceded by a radiative shock should be un- 
stable, similarly to the ionization-shock (I-S) front insta- 
bility of a D-type ionization front. By analogy to the 
I-S front instability (Giuliani (1979)), we called it here 
dissociation-shock (D-S) front instability. The dynamical be- 
havior of the D-S front instability is equivalent to the one 
of the I-S front instability ( jGarcfa- Segura & Franco (1996) 



|Whalen & Nor man (2008)), the only difference is that the tem 
perature of a PDR is around one order of magnitude lower than 
that of a HII region. This causes a lower pressure gradient be- 
tween the molecular and the dissociated gas, and consequently, 
the growth, development, and dynamics of the D-S front insta- 
bility would be less violent than that of the I-S front instability. 

At the onset of the D-S front instability, the CSM 
gas is forced to cross oblique shocks, and the gas ac- 
quires a transverse component in velocity (Vish niac (1983)| 
|Mac Low & Norman (19 93)). Thus, the gas piles up in the "val- 
leys" (by analogy), enhancing the column density, but the gas 
is cleared out from the "peaks", therefore lowering the column 
density. The radiation field notices the redistribution of the gas 
and a number of photons can now reach deeper regions in the 
directions where the "peaks" are located. Because the inertia of 
the swept-up shell decreases at the "peaks", the dissociated gas 
expands and accelerates the leading shock, steepening the an- 
gle of the oblique shocks. If the CSM has a decreasing slope in 
density (such as r~ 2 in AGB winds), the growth of this insta- 
bility becomes catastrophic. Dissociation fronts should be also 
subject to a shadowing instability similar to ionization fronts 
( [Williams (1999)] |Whalen & Norman (20083) , which probably 
triggers the development of the D-S front instability. 

We have shown here that photodissociation produces an im- 
portant fragmentation in PPNs via the D-S front instability, and 
is able to shape complicated structures, like quadrupolar and 
multipolar lobes. Then, according to Kwok's idea, photodisso- 
ciation is the mechanism able to create holes in the surround- 
ing medium of post-AGBs. Once the holes are formed, the post- 
AGB wind could be channeled through them, forming multipolar 
bubbles. 

In a future article, we will perform calculations using new 
evolutionary tracks for post-AGBs, including other inputs such 
as stellar winds, with resolution studies for the D-S front insta- 
bility, and with 3-D solutions. 
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